1. Field of the Invention
The present invention relates to a method for measuring an internal property distribution in a scattering medium and an apparatus therefor. More specifically, the present invention concerns a measuring method of internal property distribution and is applicable to equipment for obtaining internal information while moving a light injection position and a light detection position along a surface of a measured object, and an apparatus therefor.
2. Related Background Art
Optical CT (computer tomography) means a technique or apparatus of measuring an optical property distribution or a concentration distribution of an absorptive constituent in an organism, and makes use of light (signal light) detected after injection of light into a living tissue and migration therethrough.
Three types of optical CT techniques making use of rectilinearly propagating light, quasi-rectilinearly propagating light, and scattered light have been reported heretofore. Among these, the method making use of rectilinearly propagating light has extremely poor utilization efficiency of light and is thus applicable only to very small media. For example, where the near-infrared light is used, transmitted light demonstrates attenuation to about 10xe2x88x925 times against a standard living tissue having the thickness of 1 cm. In contrast with it, the method making use of scattered light utilizes all light emerging from the medium, so as to increase signal-to-noise ratios, and is thus expected to be applied to larger media. The method making use of quasi-rectilinearly propagating light stands intermediate between these two methods.
The practical optical CT needs to utilize the scattered light because of the limitations including maximum permissible incidence power to organisms, measurement sensitivity, required measurement time, and so on, but the above optical CT does not have been put into practical use yet because of the following technological problems.
The first problem is that no method has been developed for describing the behavior of light or a photon in a scattering medium with sufficient accuracy. For analyzing the behavior of a photon migrating in a scattering medium, it has been common practice heretofore to employ approximation to the transport equation, or the photon diffusion equation resulting from application of diffusion approximation to the transport theory. The diffusion approximation, however, holds only in sufficiently larger media than the mean free path length of photons therein and is thus incapable of handling relatively small media, tissues having complicated internal shapes, and media having complicated shapes. In addition, the diffusion approximation is predicated on isotropic scattering; therefore, when applied to measurement of actual living tissues having anisotropic scattering characteristics, it gives rise to unignorable errors due to the anisotropic scattering. Further, the diffusion equation does not allow us to find its solution by either of analytical or numerical techniques (such as the finite element method or the like) unless boundary conditions are preliminarily set. Namely, it is necessary to set the boundary conditions at each of the light injection and detection positions, i.e., the shape of the medium and reflection characteristics at interfaces, prior to measurement. If these conditions vary depending upon individual differences etc., computation must be redone under boundary conditions modified according to the variation. Therefore, the optical CT making use of the relation between signal light and optical properties of a scattering medium derived from the approximate expression of the transport equation or the photon diffusion equation still has a significant problem in accuracy and operability.
Besides the above, there is another method for deriving the relation between signal light and optical properties of a scattering medium by applying the perturbation theory to the approximate expression of the transport equation or to the photon diffusion equation, and for reconstructing an optical CT image using this relation. This method, however, makes how to handle its nonlinear effects (terms of the second and higher degrees) very complex. On this occasion, it is theoretically possible to perform the computation including the terms of the second and higher degrees by a computer, but operation time is enormous even with use of the presently fastest computer; therefore, practical use thereof is impossible. It is thus usual practice to ignore the terms of the second and higher degrees. Therefore, when this method is applied to reconstruction of an optical CT image of a medium containing a plurality of relatively strong absorption regions, interaction becomes unignorable between the absorption regions and a large error can be made due to it.
The second problem is that the conventional optical CT makes use of a weight function in a narrow sense, i.e., a mean path length or a phase lag equivalent thereto. For this reason, it becomes extremely complicated to handle the mean path length of detected light varying depending upon absorption coefficients. It is thus usual practice to employ approximation, but use of approximation poses a significant problem of increase in errors. Such methods making use of the weight function in a narrow sense are described, for example, in references listed below. (1) S. Arridge: SPIE Institutes for Advanced Optical Technologies, Vol. IS11, Medical Optical Tomography: Functional Imaging and Monitoring, 35-64 (1993); (2) R. L. Barbour and H. L. Graber: ibid. 87-120 (1993); (3) H. L. Graber, J. Chang, R. Aronson and R. L. Barbour: ibid. 121-143 (1993); (4) J. C. Schotland, J. C. Haselgrove and J. S. Leigh: Applied optics, 32, 448-5453 (1993); (5) Chang, R. Aronson, H. L. Graber and R. L. Barbour: Proc. SPIE, 2389, 448-464 (1995); (6) B. W. Pogue, M. S. Patterson, H. Jiang and K. D. Paulsen: Phys. Med. Biol. 40, 1709-1729 (1995); (7) S. R. Arridge: Applied Optics, 34, 7395-7409 (1995); (8) H. L. Graber, J. Chang, and R. L. Barbour: Proc. SPIE, 2570, 219-234 (1995); (9) A. Maki and H. Koizumi: OSA TOPS, Vol. 2, 299-304 (1996); (10) H. Jiang, K. D. Paulsen and Ulf L. Osterberg: J. Opt. Soc. Am. A13, 253-266 (1996); (11) S. R. Arridge and J. C. Hebden: Phys. Med. Biol. 42, 841-853 (1997); (12) S. B. Colak, D. G. Papaioannou, G. W. It Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen and N. A. A. J. van Astten: Applied Optics, 36, 180-213 (1997).
In image reconstruction of optical CT, it is most important to know which part the signal light has passed in the scattering medium, i.e., to know a path distribution of the signal light in the medium. From this viewpoint, there is also a method for deriving the path distribution of signal light by the random walk theory or the like, but the aforementioned problem is not solved yet.
As described above, the conventional optical CT techniques do not allow us to obtain a reconstructed image with sufficient accuracy and still have significant issues in terms of spatial resolution, image distortion, quantitation, measurement sensitivity, required measurement time, and so on.
In order to break through the circumstances as described above, the inventors have been conducting a series of studies with the focus on the following points. Namely, the important points for realizing optical CT are to clarify the behavior of light migrating in living tissues of strong scattering media, to clarify the relation between signal light detected and optical properties of scattering media (scattering absorbers) containing absorptive constituents, and to develop an algorithm for reconstructing an optical CT image by making use of the signal light and the relation.
Then the inventors proposed {circle around (1)} a model based on Microscopic Beer-Lambert Law (hereinafter referred to as xe2x80x9cMBLxe2x80x9d), derived {circle around (2)} analytic expressions to indicate the relation between optical properties of scattering media and signal light, and reported them in the following references, for example. (13) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 34, L79-81 (1995); (14) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 35, 4848-4851 (1996); (15) Y. Tsuchiya and Urakami: Optics Communications, 144, 269-280 (1997); (16) H. Zhang, M. Miwa, Y. Yamashita and Y. Tsuchiya: Ext. Abstr. Optics Japan ""97, 30aA08 (1997). From the principle, this MBL-based method has a significant feature of being free of influence of the medium shape, boundary conditions, and scattering and thus is also applicable to anisotropic scattering media and small media. As a consequence, the method clarified the photon migration in scattering media and made it possible to measure absorption coefficients of scattering media or concentrations of absorbers, based on this theory.
The inventors, however, found that neither of the above conventional methods was readily applicable to inhomogeneous systems, because they made use of the conventional weight function in a narrow sense, and that the measurement accuracy was not sufficient yet for the absorptive constituents in the scattering media such as organisms.
The present invention has been accomplished in view of the above conventional issues and an object of the invention is to provide a method and an apparatus that permit measurement of an absorptive constituent with higher accuracy even when applied to the scattering media of inhomogeneous systems such as organisms or the like.
As a result of extensive and intensive research in order to achieve the above object, the inventors further developed the above knowledge found by the inventor, and discovered that the above object was accomplished by {circle around (3)} applying the aforementioned MBL to the inhomogeneous systems, {circle around (4)} deriving analytic expressions to indicate the relation between optical properties of inhomogeneous scattering media and signal light, {circle around (5)} deriving a weight function according to a different definition from the conventional one from the relation between signal light and optical properties of inhomogeneous scattering media, and {circle around (6)} providing an algorithm and an apparatus for reconstructing an optical CT image by use of this weight function, whereby the inventors reached the invention.
Namely, a method for measuring an internal property distribution of a scattering medium according to the present invention is a method comprising:
a light injection step of successively injecting rays from at least one light injection position into a medium to be measured, which is a scattering medium;
a light detection step of detecting rays having passed through the interior of the medium to be measured, at a plurality of light detection positions;
a measurement value acquisition step of acquiring a measurement value of a predetermined parameter of the rays for each of combinations of the light injection position with the light detection positions, based on each ray detected;
a reference value setting step of setting a reference value of an absorption coefficient of the medium to be measured;
an estimate computation step of computing an estimate of said parameter for each of the combinations of the light injection position with the light detection positions, based on the reference value of the absorption coefficient, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole;
a weight function operation step of obtaining a weight function in each voxel of the medium to be measured, the medium being divided into a plurality of voxels, based on the Microscopic Beer-Lambert Law, using the reference value of the homogeneous absorption coefficient;
an absorption coefficient deviation computation step of computing a deviation of the absorption coefficient from the reference value of the absorption coefficient in each voxel, based on the measurement value of the parameter, the estimate of the parameter, and the weight function; and
an absorption coefficient absolute value computation step of computing an absolute value of the absorption coefficient in each voxel, based on the reference value of the absorption coefficient and the deviation of the absorption coefficient, to obtain a distribution of absolute values of absorption coefficients in the medium to be measured.
An apparatus for measuring an internal property distribution of a scattering medium according to the present invention is an apparatus comprising:
light injection means for successively injecting rays from at least one light injection position into a medium to be measured, which is a scattering medium;
light detection means for detecting rays having passed through the interior of the medium to be measured, at a plurality of light detection positions;
measurement value acquisition means for acquiring a measurement value of a predetermined parameter of the rays for each of combinations of the light injection position with the light detection positions, based on each ray detected;
reference value setting means for setting a reference value of an absorption coefficient of the medium to be measured;
estimate computation means for computing an estimate of the parameter for each of the combinations of the light injection position with the light detection positions, based on the reference value of the absorption coefficient, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole;
weight function operation means for obtaining a weight function in each voxel of the medium to be measured, the medium being divided into a plurality of voxels, based on the Microscopic Beer-Lambert Law, using the homogeneous reference value of the absorption coefficient;
absorption coefficient deviation computation means for computing a deviation of the absorption coefficient from the reference value of the absorption coefficient in each voxel, based on the measurement value of the parameter, the estimate of the parameter, and the weight function; and
absorption coefficient absolute value computation means for computing an absolute value of the absorption coefficient in each voxel, based on the reference value of the absorption coefficient and the deviation of the absorption coefficient, to obtain a distribution of absolute values of absorption coefficients in the medium to be measured.
In the method and apparatus of the present invention, in each measurement the weight function in each voxel, which is first disclosed by the present invention, is directly gained based on the MBL, and the deviation of the absorption coefficient in each voxel is computed based on the weight function, the measurement value of the predetermined parameter, and the estimate of the parameter. On this occasion, the weight function is expressed by an equation of one kind regardless of the measurement circumstances. Since in the present invention the deviation of the absorption coefficient is computed based on the appropriate weight function even with variations in the individual measurement circumstances as described above, the invention prevents errors from being caused by employment of the approximation. For this reason, even in cases where the method and apparatus of the present invention are applied to the scattering media of inhomogeneous systems such as organisms or the like, they allow the deviation of the absorption coefficient in each voxel to be computed accurately, and improve the measurement accuracy of the internal property distribution (optical CT image) obtained based on such absorption coefficient deviation.
A preferred weight function adopted in the above method and apparatus of the present invention is a function of the mean path length in each voxel and the variance of the distribution of path lengths, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole.
In this case, it is preferable that the method (or the apparatus) further comprise the mean path length acquisition step (or the mean path length acquisition means) of acquiring the mean path length in each voxel, based on the reference value of the absorption coefficient, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole, and that the aforementioned weight function operation step comprise (or the weight function operation means perform) a step of gaining the weight function, based on the following equation:
Wi=Zi(xcexcav)xe2x88x92(xcex94xcexcai/2)"sgr"iv2xe2x80x83xe2x80x83(1)
[where Wi is the weight function, Zi(xcexcav) the mean path length, xcex94xcexcai the deviation of the absorption coefficient, and "sgr"iv2 the variance of the distribution of path lengths].
Another preferred weight function adopted in the above method and apparatus of the present invention is a function of the mean path length in a predetermined time domain in each voxel and a variance of a distribution of path lengths, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole.
In this case, it is preferable that the method (or the apparatus) further comprise the mean path length acquisition step (or the mean path length acquisition means) of obtaining the mean path length in a predetermined time domain in each voxel, based on the reference value of the absorption coefficient, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole, and that the aforementioned weight function operation step comprise (or the weight function operation means perform) a step of gaining the weight function, based on the following equation:
Wi=Zi(xcexcav)xe2x88x92(xcex94xcexcai/2)"sgr"iv2xe2x80x83xe2x80x83(1 )
[where Wi is the weight function, Zi(xcexcav) the mean path length, xcex94xcexcai the deviation of the absorption coefficient, and "sgr"iv2 the variance of the distribution of path lengths].
Further, still another preferred weight function adopted in the above method and apparatus of the present invention is a function of a group delay in each voxel and a variance of a distribution thereof, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole.
In this case, it is preferable that the method (or the apparatus) further comprise a group delay acquisition step (or group delay acquisition means) of acquiring a group delay in each voxel, based on the reference value of the absorption coefficient, on the assumption that the medium to be measured has the homogeneous reference value of the absorption coefficient as a whole, and that the aforementioned weight function operation step comprise (or the weight function operation means perform) a step of gaining the weight function, based on the following equation:                                                         "AutoLeftMatch"                                                W                  i                                =                                  c                  ⁢                                      xe2x80x83                                    ⁢                                                            ∂                                              φ                        i                                                                                    ∂                      ω                                                                                  "RightBracketingBar"                        "AutoRightMatch"                                μ            av                          -                                            Δ              ⁢                              xe2x80x83                            ⁢                              μ                ai                                      2                    ⁢                      xe2x80x83                    ⁢                      σ            f            2                                              (        2        )            
[where Wi is the weight function,             "AutoLeftMatch"              c        ⁢                  xe2x80x83                ⁢                              ∂                          φ              i                                            ∂            ω                              "RightBracketingBar"        "AutoRightMatch"        μ    av  
is c (speed of light in the medium) times the group delay, xcex94xcexcai the deviation of the absorption coefficient, and "sgr"f2 the variance of a distribution].
With adoption of the preferred weight functions in the method and apparatus of the present invention as described above, because the weight functions are computed in consideration of not only the mean path length or the group delay varying depending upon the absorption coefficient, but also the variance of distribution thereof, there is a tendency to attain further improvement in the measurement accuracy of the internal property distribution obtained by use of such weight functions.
The method and apparatus of the present invention may further comprise a concentration computation step (or concentration computation means) of computing a concentration of an absorptive constituent in each voxel by using the absolute value of the absorption coefficient, and thus obtaining a concentration distribution of the absorption constituent in the medium to be measured. Since concentrations of the absorptive constituent can be obtained based on the absorption coefficient deviations determined accurately as described above by such method and apparatus, the concentration distribution can be obtained with high accuracy.
In the cases where the method and apparatus of the present invention are applied to a measured medium containing at least two absorptive constituents, it is preferable that the rays injected into the measured medium in the light injection step (light injection means) have at least two wavelengths at which said absorptive constituents demonstrate their respective absorption coefficients different from each other. In this case, each of the rays having the at least two wavelengths is detected in the light detection step (or by the light detection means), the measurement value is acquired as to each of the rays having the at least two wavelengths in the measurement value acquisition step (or by the measurement value acquisition means), the reference value is set as to each of the rays having the at least two wavelengths in the reference value setting step (or by the reference value setting means), the estimate is computed as to each of the rays having the at least two wavelengths in the estimate computation step (or by the estimate computation means), the weight function is gained as to each of the rays having the at least two wavelengths in the weight function operation step (or by the weight function operation means), the deviation of the absorption coefficient is computed as to each of the rays having the at least two wavelengths in the absorption coefficient deviation computation step (or by the absorption coefficient deviation computation means), the absolute value of the absorption coefficient is computed as to each of the rays having the at least two wavelengths in the absorption coefficient absolute value computation step (or by the absorption coefficient absolute value computation means), and concentrations of each absorptive constituent are computed as to each of the rays having the at least two wavelengths in the concentration computation step (or by the concentration computation means), thereby obtaining a concentration distribution of each of the absorption constituents in the measured medium.
The method and apparatus of the present invention described above may further comprise an image display step (or image display means) of displaying an optical CT image to indicate the distribution in the measured medium, based on the aforementioned distribution obtained. The method and apparatus of the present invention in this structure can display the optical CT image with high accuracy by imaging of the internal property distribution obtained with accuracy as described above.